USGS has libraries to retreive data from their system.

install libraries - do only one time when first runnign script

# https://github.com/USGS-R/dataRetrieval

# general R packages to install if you dont have them
# remove the "#" if you need to run the installations

# install.packages("devtools") # essential in installing other thigs
# install.packages("tidyverse") # installs a lot of things and ggplot
# install.packages("scales") # allows great scale formatting on ggplot
# install.packages("lubridate") # makes working with dates easier
# install.packages("janitor") # clean names of columns and other things
# install.packages("readxl") # read in excel files
# install.packages("plotly") # interactive plots
# devtools::install_github("thomasp85/patchwork") # multiple plots

# specific to this module
# install.packages("dataRetrieval") # USGS Data Retreiveal Method

Libraires

library(skimr)

Attaching package: ‘skimr’

The following object is masked from ‘package:knitr’:

    kable

The following object is masked from ‘package:stats’:

    filter

Read in file for …


# The site is Neversink and the site number is 01435000
# USGS 01435000 NEVERSINK RIVER NEAR CLARYVILLE NY

siteNo_monthly <-     "01435000"
pCode_monthly <-      "00060"
start.date_monthly <- "1937-11"
end.date_monthly <-   "2019-08"

neversink_monthly.df <-   readNWISstat(siteNumbers = siteNo_monthly,
                        parameterCd = pCode_monthly,
                        startDate = start.date_monthly,
                        endDate = end.date_monthly,
                        statReportType="Monthly",
                        statType="mean")
Please be aware the NWIS data service feeding this function is in BETA.

          Data formatting could be changed at any time, and is not guaranteed
  
# rename the columns
neversink_monthly.df <- renameNWISColumns(neversink_monthly.df)

names(neversink_monthly.df)
[1] "agency_cd"    "site_no"      "parameter_cd" "ts_id"        "loc_web_ds"  
[6] "year_nu"      "month_nu"     "mean_va"     

we need to make a date column to plot with

neversink_monthly.df <- neversink_monthly.df %>%
  mutate(
    date = ymd(paste(year_nu, month_nu, "01", sep="-"))
  )

Plot all data

all.plot <- ggplot(neversink_monthly.df, aes(date, mean_va)) +
  geom_point() +
  geom_line()
all.plot

Plot February data

feb.plot <- neversink_monthly.df %>%
  filter(month_nu == 2) %>%
  ggplot(aes(date, mean_va)) +
  geom_point() +
  geom_line() +
  labs(y="Mean Monthly Discharge m^3/sec", x="Date")
feb.plot

NA

We can make it interactive

ggplotly(feb.plot)

August data

aug.plot <- neversink_monthly.df %>%
  filter(month_nu == 8) %>%
  ggplot(aes(date, mean_va)) +
  geom_point() +
  geom_line() +
  labs(y="Mean Monthly Discharge m^3/sec", x="Date")
aug.plot

Compare the plots

feb_aug.plot <- feb.plot +
                aug.plot +
                plot_layout(ncol =1)
feb_aug.plot

add trend lines as straight lines

feb_aug.plot <- feb.plot + geom_smooth(method="lm") +
                aug.plot + geom_smooth(method="lm") +
                plot_layout(ncol =1)
feb_aug.plot

How to see regression statistics

summary(feb.model)

Call:
lm(formula = mean_va ~ date, data = .)

Residuals:
   Min     1Q Median     3Q    Max 
-184.1 -114.4  -36.8   76.0  697.1 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 1.944e+02  5.148e+00   37.76   <2e-16 ***
date        8.601e-04  5.660e-04    1.52    0.129    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 147.5 on 943 degrees of freedom
Multiple R-squared:  0.002443,  Adjusted R-squared:  0.001385 
F-statistic: 2.309 on 1 and 943 DF,  p-value: 0.1289
summary(aug.model)

Call:
lm(formula = mean_va ~ date, data = .)

Residuals:
   Min     1Q Median     3Q    Max 
-184.1 -114.4  -36.8   76.0  697.1 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 1.944e+02  5.148e+00   37.76   <2e-16 ***
date        8.601e-04  5.660e-04    1.52    0.129    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 147.5 on 943 degrees of freedom
Multiple R-squared:  0.002443,  Adjusted R-squared:  0.001385 
F-statistic: 2.309 on 1 and 943 DF,  p-value: 0.1289

Try another site

Use this site to find a new location https://maps.waterdata.usgs.gov/mapper/index.html

# The site is Neversink and the site number is 01435000
# USGS SAGAVANIRKTOK R NR PUMP STA 3 AK

siteNo_monthly <-     "15908000"
pCode_monthly <-      "00060"
start.date_monthly <- "1900-01"
end.date_monthly <-   "2019-08"

sag_monthly.df <-   readNWISstat(siteNumbers = siteNo_monthly,
                        parameterCd = pCode_monthly,
                        startDate = start.date_monthly,
                        endDate = end.date_monthly,
                        statReportType="Monthly",
                        statType="mean")
Please be aware the NWIS data service feeding this function is in BETA.

          Data formatting could be changed at any time, and is not guaranteed
  
# rename the columns
sag_monthly.df <- renameNWISColumns(sag_monthly.df)

names(sag_monthly.df)
[1] "agency_cd"    "site_no"      "parameter_cd" "ts_id"        "loc_web_ds"  
[6] "year_nu"      "month_nu"     "mean_va"     

we need to make a date column to plot with

sag_monthly.df <- sag_monthly.df %>%
  mutate(
    date = ymd(paste(year_nu, month_nu, "01", sep="-"))
  )

Plot February data

feb_new.plot <- sag_monthly.df %>%
  filter(month_nu == 2) %>%
  ggplot(aes(date, mean_va)) +
  geom_point() +
  geom_line() +
  labs(y="Mean Monthly Discharge m^3/sec", x="Date")
feb_new.plot

NA

August data

aug_new.plot <- sag_monthly.df %>%
  filter(month_nu == 8) %>%
  ggplot(aes(date, mean_va)) +
  geom_point() +
  geom_line() +
  labs(y="Mean Monthly Discharge m^3/sec", x="Date")
aug_new.plot

Compare the plots

feb_aug_new.plot <- feb_new.plot +
                aug_new.plot +
                plot_layout(ncol =1)
feb_aug_new.plot

Never sink Peak Flows

siteNo_peak <-     "01435000"
pCode_peak <-      "00060"
start.date_peak <- "1850-01-01"
end.date_peak <-   "2019-08-02"

neversink_peak.df <-   readNWISpeak(siteNumbers = siteNo_peak)

# rename the columns
neversink_peak.df <- renameNWISColumns(neversink_peak.df)

names(neversink_peak.df)

Plot the data

peak_flow.plot <- neversink_peak.df  %>%
  ggplot(aes(peak_dt , peak_va)) + 
  geom_point() +
  geom_line()
ggplotly(peak_flow.plot)

mean and std flows

skim(neversink_peak.df)
Skim summary statistics
 n obs: 79 
 n variables: 13 

── Variable type:character ────────────────────────────────────────────────────────────────────────────────────────────────────
      variable missing complete  n min max empty n_unique
 ag_gage_ht_cd      76        3 79   1   3     0        2
         ag_tm      79        0 79  NA  NA     0        0
     agency_cd       0       79 79   4   4     0        1
    gage_ht_cd      37       42 79   1   3     0        4
       peak_cd      73        6 79   1   1     0        3
       peak_tm      79        0 79  NA  NA     0        0
       site_no       0       79 79   8   8     0        1

── Variable type:Date ─────────────────────────────────────────────────────────────────────────────────────────────────────────
 variable missing complete  n        min        max     median n_unique
    ag_dt      76        3 79 1948-02-14 1976-01-27 1960-04-04        3
  peak_dt       0       79 79 1938-07-22 2017-02-25 1977-10-01       79

── Variable type:numeric ──────────────────────────────────────────────────────────────────────────────────────────────────────
     variable missing complete  n    mean      sd      p0     p25     p50     p75     p100     hist
   ag_gage_ht      76        3 79    6.33    0.59    5.99    6       6       6.5      7.01 ▇▁▁▁▁▁▁▃
      gage_ht       0       79 79    8.54    2.99    3.73    6.02    8.51   10.95    14.64 ▆▇▅▂▆▆▅▂
      peak_va       0       79 79 7189.24 4497.07 1400    4215    6080    9090    23400    ▆▇▅▃▁▁▁▁
 year_last_pk      78        1 79 1938      NA    1938    1938    1938    1938     1938    ▁▁▁▇▁▁▁▁

calculate flood probability

this is one way to do it but is not the best


neversink_floodfreq.df <- neversink_peak.df %>%
  select(agency_cd, site_no, peak_dt, peak_va, gage_ht) %>%
  arrange(desc(peak_va)) %>% 
  mutate(
    year =  year(peak_dt)
  ) %>%
  group_by(year) %>%
  filter(peak_va == max(peak_va)) %>%
  mutate(rank = dense_rank(desc(peak_va)))

neversink_floodfreq.df <- neversink_floodfreq.df %>%  
  ungroup() %>% 
  mutate(total_n = n()) %>% 
  rowwise() %>%
  mutate(probability = rank/(total_n+1))

LS0tCnRpdGxlOiAiRWRkaWUgbW9kdWxlIHN0cmVhbSBkaXNjaGFyZ2UgZnJvbSBoYW5kb3V0IgpvdXRwdXQ6IGh0bWxfbm90ZWJvb2sKZWRpdG9yX29wdGlvbnM6IAogIGNodW5rX291dHB1dF90eXBlOiBpbmxpbmUKLS0tCgojIFVTR1MgaGFzIGxpYnJhcmllcyB0byByZXRyZWl2ZSBkYXRhIGZyb20gdGhlaXIgc3lzdGVtLiAgICAKCiMgaW5zdGFsbCBsaWJyYXJpZXMgLSBkbyBvbmx5IG9uZSB0aW1lIHdoZW4gZmlyc3QgcnVubmlnbiBzY3JpcHQgICAgCmBgYHtyIGluc3RhbGwgcGFja2FnZXMsIG1lc3NhZ2U9VFJVRSwgd2FybmluZz1UUlVFfQojIGh0dHBzOi8vZ2l0aHViLmNvbS9VU0dTLVIvZGF0YVJldHJpZXZhbAoKIyBnZW5lcmFsIFIgcGFja2FnZXMgdG8gaW5zdGFsbCBpZiB5b3UgZG9udCBoYXZlIHRoZW0KIyByZW1vdmUgdGhlICIjIiBpZiB5b3UgbmVlZCB0byBydW4gdGhlIGluc3RhbGxhdGlvbnMKCiMgaW5zdGFsbC5wYWNrYWdlcygiZGV2dG9vbHMiKSAjIGVzc2VudGlhbCBpbiBpbnN0YWxsaW5nIG90aGVyIHRoaWdzCiMgaW5zdGFsbC5wYWNrYWdlcygidGlkeXZlcnNlIikgIyBpbnN0YWxscyBhIGxvdCBvZiB0aGluZ3MgYW5kIGdncGxvdAojIGluc3RhbGwucGFja2FnZXMoInNjYWxlcyIpICMgYWxsb3dzIGdyZWF0IHNjYWxlIGZvcm1hdHRpbmcgb24gZ2dwbG90CiMgaW5zdGFsbC5wYWNrYWdlcygibHVicmlkYXRlIikgIyBtYWtlcyB3b3JraW5nIHdpdGggZGF0ZXMgZWFzaWVyCiMgaW5zdGFsbC5wYWNrYWdlcygiamFuaXRvciIpICMgY2xlYW4gbmFtZXMgb2YgY29sdW1ucyBhbmQgb3RoZXIgdGhpbmdzCiMgaW5zdGFsbC5wYWNrYWdlcygicmVhZHhsIikgIyByZWFkIGluIGV4Y2VsIGZpbGVzCiMgaW5zdGFsbC5wYWNrYWdlcygicGxvdGx5IikgIyBpbnRlcmFjdGl2ZSBwbG90cwojIGluc3RhbGwucGFja2FnZXMoc2tpbXIpCiMgZGV2dG9vbHM6Omluc3RhbGxfZ2l0aHViKCJ0aG9tYXNwODUvcGF0Y2h3b3JrIikgIyBtdWx0aXBsZSBwbG90cwoKIyBzcGVjaWZpYyB0byB0aGlzIG1vZHVsZQojIGluc3RhbGwucGFja2FnZXMoImRhdGFSZXRyaWV2YWwiKSAjIFVTR1MgRGF0YSBSZXRyZWl2ZWFsIE1ldGhvZApgYGAKCgpMaWJyYWlyZXMKYGBge3IgbGFvZCBsaWJyYXJpZXN9CiMgbG9hZCB0aGVzZSBldmVyeSB0aW1lICAKbGlicmFyeSh0aWR5dmVyc2UpCmxpYnJhcnkoc2NhbGVzKQpsaWJyYXJ5KGx1YnJpZGF0ZSkKbGlicmFyeShqYW5pdG9yKQpsaWJyYXJ5KHJlYWR4bCkKbGlicmFyeShza2ltcikKbGlicmFyeShwbG90bHkpCmxpYnJhcnkocGF0Y2h3b3JrKQoKIyBsb2FkIHRoZSBVU0dTIHBhY2thZ2UKbGlicmFyeShkYXRhUmV0cmlldmFsKQpgYGAKCgoKUmVhZCBpbiBmaWxlIGZvciAuLi4KYGBge3IgcmVhZCBpbiBmaWxlfQoKIyBUaGUgc2l0ZSBpcyBOZXZlcnNpbmsgYW5kIHRoZSBzaXRlIG51bWJlciBpcyAwMTQzNTAwMAojIFVTR1MgMDE0MzUwMDAgTkVWRVJTSU5LIFJJVkVSIE5FQVIgQ0xBUllWSUxMRSBOWQoKc2l0ZU5vX21vbnRobHkgPC0gICAgICIwMTQzNTAwMCIKcENvZGVfbW9udGhseSA8LSAgICAgICIwMDA2MCIKc3RhcnQuZGF0ZV9tb250aGx5IDwtICIxOTM3LTExIgplbmQuZGF0ZV9tb250aGx5IDwtICAgIjIwMTktMDgiCgpuZXZlcnNpbmtfbW9udGhseS5kZiA8LSAgIHJlYWROV0lTc3RhdChzaXRlTnVtYmVycyA9IHNpdGVOb19tb250aGx5LAogICAgICAgICAgICAgICAgICAgICAgICBwYXJhbWV0ZXJDZCA9IHBDb2RlX21vbnRobHksCiAgICAgICAgICAgICAgICAgICAgICAgIHN0YXJ0RGF0ZSA9IHN0YXJ0LmRhdGVfbW9udGhseSwKICAgICAgICAgICAgICAgICAgICAgICAgZW5kRGF0ZSA9IGVuZC5kYXRlX21vbnRobHksCiAgICAgICAgICAgICAgICAgICAgICAgIHN0YXRSZXBvcnRUeXBlPSJNb250aGx5IiwKICAgICAgICAgICAgICAgICAgICAgICAgc3RhdFR5cGU9Im1lYW4iKQogIAojIHJlbmFtZSB0aGUgY29sdW1ucwpuZXZlcnNpbmtfbW9udGhseS5kZiA8LSByZW5hbWVOV0lTQ29sdW1ucyhuZXZlcnNpbmtfbW9udGhseS5kZikKCm5hbWVzKG5ldmVyc2lua19tb250aGx5LmRmKQpgYGAKCiMgd2UgbmVlZCB0byBtYWtlIGEgZGF0ZSBjb2x1bW4gdG8gcGxvdCB3aXRoCmBgYHtyfQpuZXZlcnNpbmtfbW9udGhseS5kZiA8LSBuZXZlcnNpbmtfbW9udGhseS5kZiAlPiUKICBtdXRhdGUoCiAgICBkYXRlID0geW1kKHBhc3RlKHllYXJfbnUsIG1vbnRoX251LCAiMDEiLCBzZXA9Ii0iKSkKICApCgpgYGAKCgoKIyBQbG90IGFsbCBkYXRhCmBgYHtyfQphbGwucGxvdCA8LSBnZ3Bsb3QobmV2ZXJzaW5rX21vbnRobHkuZGYsIGFlcyhkYXRlLCBtZWFuX3ZhKSkgKwogIGdlb21fcG9pbnQoKSArCiAgZ2VvbV9saW5lKCkKYWxsLnBsb3QKYGBgCgoKCiMgUGxvdCBGZWJydWFyeSBkYXRhIAoKYGBge3J9CmZlYi5wbG90IDwtIG5ldmVyc2lua19tb250aGx5LmRmICU+JQogIGZpbHRlcihtb250aF9udSA9PSAyKSAlPiUKICBnZ3Bsb3QoYWVzKGRhdGUsIG1lYW5fdmEpKSArCiAgZ2VvbV9wb2ludCgpICsKICBnZW9tX2xpbmUoKSArCiAgbGFicyh5PSJNZWFuIE1vbnRobHkgRGlzY2hhcmdlIG1eMy9zZWMiLCB4PSJEYXRlIikKZmViLnBsb3QKICAKYGBgCgpXZSBjYW4gbWFrZSBpdCBpbnRlcmFjdGl2ZQpgYGB7cn0KZ2dwbG90bHkoZmViLnBsb3QpCmBgYAoKCkF1Z3VzdCBkYXRhCmBgYHtyfQphdWcucGxvdCA8LSBuZXZlcnNpbmtfbW9udGhseS5kZiAlPiUKICBmaWx0ZXIobW9udGhfbnUgPT0gOCkgJT4lCiAgZ2dwbG90KGFlcyhkYXRlLCBtZWFuX3ZhKSkgKwogIGdlb21fcG9pbnQoKSArCiAgZ2VvbV9saW5lKCkgKwogIGxhYnMoeT0iTWVhbiBNb250aGx5IERpc2NoYXJnZSBtXjMvc2VjIiwgeD0iRGF0ZSIpCmF1Zy5wbG90CmBgYAoKIyBDb21wYXJlIHRoZSBwbG90cwpgYGB7cn0KZmViX2F1Zy5wbG90IDwtIGZlYi5wbG90ICsKICAgICAgICAgICAgICAgIGF1Zy5wbG90ICsKICAgICAgICAgICAgICAgIHBsb3RfbGF5b3V0KG5jb2wgPTEpCmZlYl9hdWcucGxvdApgYGAKCiMgYWRkIHRyZW5kIGxpbmVzIGFzIHN0cmFpZ2h0IGxpbmVzCmBgYHtyfQpmZWJfYXVnLnBsb3QgPC0gZmViLnBsb3QgKyBnZW9tX3Ntb290aChtZXRob2Q9ImxtIikgKwogICAgICAgICAgICAgICAgYXVnLnBsb3QgKyBnZW9tX3Ntb290aChtZXRob2Q9ImxtIikgKwogICAgICAgICAgICAgICAgcGxvdF9sYXlvdXQobmNvbCA9MSkKZmViX2F1Zy5wbG90CmBgYAoKIyBIb3cgdG8gc2VlIHJlZ3Jlc3Npb24gc3RhdGlzdGljcwpgYGB7cn0KZmViLm1vZGVsIDwtIG5ldmVyc2lua19tb250aGx5LmRmICU+JSAKICBmaWx0ZXIobW9udGhfbnUgPT0gMikgJT4lCiAgbG0obWVhbl92YSB+IGRhdGUsIGRhdGE9LikKc3VtbWFyeShmZWIubW9kZWwpCgpgYGAKCgpgYGB7cn0KYXVnLm1vZGVsIDwtIG5ldmVyc2lua19tb250aGx5LmRmICU+JSAKICBmaWx0ZXIobW9udGhfbnUgPT0gOCkgJT4lCiAgbG0obWVhbl92YSB+IGRhdGUsIGRhdGE9LikKc3VtbWFyeShhdWcubW9kZWwpCgpgYGAKCgoKIyBUcnkgYW5vdGhlciBzaXRlClVzZSB0aGlzIHNpdGUgdG8gZmluZCBhIG5ldyBsb2NhdGlvbgpodHRwczovL21hcHMud2F0ZXJkYXRhLnVzZ3MuZ292L21hcHBlci9pbmRleC5odG1sCgpgYGB7cn0KIyBUaGUgc2l0ZSBpcyBOZXZlcnNpbmsgYW5kIHRoZSBzaXRlIG51bWJlciBpcyAwMTQzNTAwMAojIFVTR1MgU0FHQVZBTklSS1RPSyBSIE5SIFBVTVAgU1RBIDMgQUsKCnNpdGVOb19tb250aGx5IDwtICAgICAiMTU5MDgwMDAiCnBDb2RlX21vbnRobHkgPC0gICAgICAiMDAwNjAiCnN0YXJ0LmRhdGVfbW9udGhseSA8LSAiMTk4Mi0wMSIKZW5kLmRhdGVfbW9udGhseSA8LSAgICIyMDE5LTA4IgoKc2FnX21vbnRobHkuZGYgPC0gICByZWFkTldJU3N0YXQoc2l0ZU51bWJlcnMgPSBzaXRlTm9fbW9udGhseSwKICAgICAgICAgICAgICAgICAgICAgICAgcGFyYW1ldGVyQ2QgPSBwQ29kZV9tb250aGx5LAogICAgICAgICAgICAgICAgICAgICAgICBzdGFydERhdGUgPSBzdGFydC5kYXRlX21vbnRobHksCiAgICAgICAgICAgICAgICAgICAgICAgIGVuZERhdGUgPSBlbmQuZGF0ZV9tb250aGx5LAogICAgICAgICAgICAgICAgICAgICAgICBzdGF0UmVwb3J0VHlwZT0iTW9udGhseSIsCiAgICAgICAgICAgICAgICAgICAgICAgIHN0YXRUeXBlPSJtZWFuIikKICAKIyByZW5hbWUgdGhlIGNvbHVtbnMKc2FnX21vbnRobHkuZGYgPC0gcmVuYW1lTldJU0NvbHVtbnMoc2FnX21vbnRobHkuZGYpCgpuYW1lcyhzYWdfbW9udGhseS5kZikKYGBgCgojIHdlIG5lZWQgdG8gbWFrZSBhIGRhdGUgY29sdW1uIHRvIHBsb3Qgd2l0aApgYGB7cn0Kc2FnX21vbnRobHkuZGYgPC0gc2FnX21vbnRobHkuZGYgJT4lCiAgbXV0YXRlKAogICAgZGF0ZSA9IHltZChwYXN0ZSh5ZWFyX251LCBtb250aF9udSwgIjAxIiwgc2VwPSItIikpCiAgKQoKYGBgCgoKIyBQbG90IEZlYnJ1YXJ5IGRhdGEgCgpgYGB7cn0KZmViX25ldy5wbG90IDwtIHNhZ19tb250aGx5LmRmICU+JQogIGZpbHRlcihtb250aF9udSA9PSAyKSAlPiUKICBnZ3Bsb3QoYWVzKGRhdGUsIG1lYW5fdmEpKSArCiAgZ2VvbV9wb2ludCgpICsKICBnZW9tX2xpbmUoKSArCiAgbGFicyh5PSJNZWFuIE1vbnRobHkgRGlzY2hhcmdlIG1eMy9zZWMiLCB4PSJEYXRlIikKZmViX25ldy5wbG90CiAgCmBgYAoKCgpBdWd1c3QgZGF0YQpgYGB7cn0KYXVnX25ldy5wbG90IDwtIHNhZ19tb250aGx5LmRmICU+JQogIGZpbHRlcihtb250aF9udSA9PSA4KSAlPiUKICBnZ3Bsb3QoYWVzKGRhdGUsIG1lYW5fdmEpKSArCiAgZ2VvbV9wb2ludCgpICsKICBnZW9tX2xpbmUoKSArCiAgbGFicyh5PSJNZWFuIE1vbnRobHkgRGlzY2hhcmdlIG1eMy9zZWMiLCB4PSJEYXRlIikKYXVnX25ldy5wbG90CmBgYAoKCiMgQ29tcGFyZSB0aGUgcGxvdHMKYGBge3J9CmZlYl9hdWdfbmV3LnBsb3QgPC0gZmViX25ldy5wbG90ICsKICAgICAgICAgICAgICAgIGF1Z19uZXcucGxvdCArCiAgICAgICAgICAgICAgICBwbG90X2xheW91dChuY29sID0xKQpmZWJfYXVnX25ldy5wbG90CmBgYAoKCgoKIyBOZXZlciBzaW5rIFBlYWsgRmxvd3MKYGBge3J9CnNpdGVOb19wZWFrIDwtICAgICAiMDE0MzUwMDAiCnBDb2RlX3BlYWsgPC0gICAgICAiMDAwNjAiCnN0YXJ0LmRhdGVfcGVhayA8LSAiMTg1MC0wMS0wMSIKZW5kLmRhdGVfcGVhayA8LSAgICIyMDE5LTA4LTAyIgoKbmV2ZXJzaW5rX3BlYWsuZGYgPC0gICByZWFkTldJU3BlYWsoc2l0ZU51bWJlcnMgPSBzaXRlTm9fcGVhaykKCiMgcmVuYW1lIHRoZSBjb2x1bW5zCm5ldmVyc2lua19wZWFrLmRmIDwtIHJlbmFtZU5XSVNDb2x1bW5zKG5ldmVyc2lua19wZWFrLmRmKQoKbmFtZXMobmV2ZXJzaW5rX3BlYWsuZGYpCgpgYGAKCgoKIyBQbG90IHRoZSBkYXRhCgpgYGB7cn0KcGVha19mbG93LnBsb3QgPC0gbmV2ZXJzaW5rX3BlYWsuZGYgICU+JQogIGdncGxvdChhZXMocGVha19kdCAsIHBlYWtfdmEpKSArIAogIGdlb21fcG9pbnQoKSArCiAgZ2VvbV9saW5lKCkKZ2dwbG90bHkocGVha19mbG93LnBsb3QpCmBgYAoKIyBtZWFuIGFuZCBzdGQgZmxvd3MKYGBge3J9CnNraW0obmV2ZXJzaW5rX3BlYWsuZGYpCmBgYAoKIyBjYWxjdWxhdGUgZmxvb2QgcHJvYmFiaWxpdHkKdGhpcyBpcyBvbmUgd2F5IHRvIGRvIGl0IGJ1dCBpcyBub3QgdGhlIGJlc3QKYGBge3J9CgpuZXZlcnNpbmtfZmxvb2RmcmVxLmRmIDwtIG5ldmVyc2lua19wZWFrLmRmICU+JQogIHNlbGVjdChhZ2VuY3lfY2QsIHNpdGVfbm8sIHBlYWtfZHQsIHBlYWtfdmEsIGdhZ2VfaHQpICU+JQogIGFycmFuZ2UoZGVzYyhwZWFrX3ZhKSkgJT4lIAogIG11dGF0ZSgKICAgIHllYXIgPSAgeWVhcihwZWFrX2R0KQogICkgJT4lCiAgZ3JvdXBfYnkoeWVhcikgJT4lCiAgZmlsdGVyKHBlYWtfdmEgPT0gbWF4KHBlYWtfdmEpKSAlPiUKICB1bmdyb3VwKCkgJT4lCiAgbXV0YXRlKHJhbmsgPSBkZW5zZV9yYW5rKGRlc2MocGVha192YSkpKQoKYGBgCgoKYGBge3IgY2FsY3VsYXRlIHByb2JhYmlsaXR5fQoKbmV2ZXJzaW5rX2Zsb29kZnJlcS5kZiA8LSBuZXZlcnNpbmtfZmxvb2RmcmVxLmRmICU+JSAgCiAgdW5ncm91cCgpICU+JSAKICBtdXRhdGUodG90YWxfbiA9IG4oKSkgJT4lIAogIHJvd3dpc2UoKSAlPiUKICBtdXRhdGUocHJvYmFiaWxpdHkgPSByYW5rLyh0b3RhbF9uKzEpKQoKYGBgCgoKYGBge3IgY2FsY3VsYXRlIGxvbmcgdGVybSBwcm9iYWJpbGl0eX0KbmV2ZXJzaW5rX2Zsb29kZnJlcS5kZiA8LSBuZXZlcnNpbmtfZmxvb2RmcmVxLmRmICU+JSAgCiAgbXV0YXRlKAogICAgcmVjdXJyZW5jZV9pbnRldmFsID0gMS9wcm9iYWJpbGl0eQogICkKCgpgYGAKCgoKYGBge3J9Cm5ldmVyc2lua19mbG9vZGZyZXEuZGYgJT4lCiAgZ2dwbG90KGFlcyhsb2cocmVjdXJyZW5jZV9pbnRldmFsKSwgbG9nKHBlYWtfdmEpKSkgKyAKICBnZW9tX3BvaW50KCkKCgpgYGAKCg==